module topo_mod

  use fiona
  use flogger
  use spark_cubed_sphere
  use params_mod
  use math_mod
  use recon
  use static_mod
  use hope_recon_mod
  use ghost_interp_mod
  use physics_const_mod

  implicit none

  private
  
  public topo_init
  public topo_deriv_init
  
  integer topo_nlon
  integer topo_nlat
  real(r_kind) topo_dlon, topo_dlat
  real(r_kind), allocatable, dimension(:  ) :: topo_lon
  real(r_kind), allocatable, dimension(:  ) :: topo_lat
  integer(2)  , allocatable, dimension(:,:) :: topo_zs
  integer(2)  , allocatable, dimension(:,:) :: topo_landfrac

contains
  subroutine topo_init(static)
    type(static_type), intent(inout) :: static
    
    call topo_read
    call topo_regrid(static)
    call topo_final
    call topo_recon(static)
    call ghost_interp_topo(static)
    call topo_unify_edges(static)
  end subroutine topo_init
  
  subroutine topo_deriv_init(static)
    type(static_type), intent(inout) :: static
    
    call topo_recon_deriv(static)
    call ghost_interp_topo_deriv(static)
    call topo_unify_edges_deriv(static)
  end subroutine topo_deriv_init

  subroutine topo_read()

    real(r_kind), allocatable :: tmp(:)
    integer j

    call topo_final

    if (proc%is_root()) call log_notice('Use ' // trim(topo_file) // ' as topography.')

    select case (topo_type)
    case ('etopo1')
      if (planet /= 'earth') call log_error('Topography file ' // trim(topo_file) // ' is used for the Earth!')
      call fiona_init()
      call fiona_open_dataset('topo', file_path=trim(topo_file))
      call fiona_get_dim('topo', 'x', size=topo_nlon)
      call fiona_get_dim('topo', 'y', size=topo_nlat)

      allocate(topo_lon(topo_nlon))
      allocate(topo_lat(topo_nlat))
      allocate(topo_zs (topo_nlon,topo_nlat))

      call fiona_start_input('topo')
      call fiona_input('topo', 'x', topo_lon)
      call fiona_input('topo', 'y', topo_lat)
      call fiona_input('topo', 'z', topo_zs )
      call fiona_end_input('topo')

      topo_lon = topo_lon * D2R
      topo_lat = topo_lat * D2R

      ! Change longitude to [0,360].
      allocate(tmp(int(topo_nlon/2)))
      tmp = topo_lon(int(topo_nlon/2)+1:)
      topo_lon(int(topo_nlon/2)+1:) = topo_lon(1:int(topo_nlon/2)) + pi2
      topo_lon(1:int(topo_nlon/2)) = tmp
      do j = 1, topo_nlat
        tmp = topo_zs(int(topo_nlon/2)+1:,j)
        topo_zs(int(topo_nlon/2)+1:,j) = topo_zs(1:int(topo_nlon/2),j)
        topo_zs(1:int(topo_nlon/2),j) = tmp
      end do
      deallocate(tmp)
    case ('gmted2010')
      if (planet /= 'earth') call log_error('Topography file ' // trim(topo_file) // ' is used for the Earth!')
      call fiona_init()
      call fiona_open_dataset('topo', file_path=trim(topo_file))
      call fiona_get_dim('topo', 'lon', size=topo_nlon)
      call fiona_get_dim('topo', 'lat', size=topo_nlat)

      allocate(topo_lon      (topo_nlon          ))
      allocate(topo_lat      (          topo_nlat))
      allocate(topo_zs       (topo_nlon,topo_nlat))
      allocate(topo_landfrac (topo_nlon,topo_nlat))

      call fiona_start_input('topo')
      call fiona_input('topo', 'lon'      , topo_lon      )
      call fiona_input('topo', 'lat'      , topo_lat      )
      call fiona_input('topo', 'htopo'    , topo_zs       )
      call fiona_input('topo', 'landfract', topo_landfrac )
      call fiona_end_input('topo')

      topo_lon = topo_lon * D2R
      topo_lat = topo_lat * D2R
    case ('mola32')
      if (planet /= 'mars') call log_error('Topography file ' // trim(topo_file) // ' is used for the Mars!')
      call fiona_init()
      call fiona_open_dataset('topo', file_path=trim(topo_file))
      call fiona_get_dim('topo', 'longitude', size=topo_nlon)
      call fiona_get_dim('topo', 'latitude' , size=topo_nlat)

      allocate(topo_lon(topo_nlon))
      allocate(topo_lat(topo_nlat))
      allocate(topo_zs (topo_nlon,topo_nlat))

      call fiona_start_input('topo')
      call fiona_input('topo', 'longitude', topo_lon)
      call fiona_input('topo', 'latitude' , topo_lat)
      call fiona_input('topo', 'alt'      , topo_zs )
      call fiona_end_input('topo')

      ! Reverse latitude dimension from the South Pole to the North Pole.
      topo_lat = topo_lat(topo_nlat:1:-1)
      topo_zs  = topo_zs (:,topo_nlat:1:-1)
    case default
      call log_error('Unknown topo_type "' // trim(topo_type) // '"!', pid=proc%id)
    end select

  end subroutine topo_read

  subroutine topo_regrid(static)
    type(static_type), intent(inout) :: static

    integer idom, i, j, ig, jg, n
    integer pvs, pve, ims, ime, jms, jme, its, ite, jts, jte
    integer irs, ire, jrs, jre
    real(r_kind) min_lon, max_lon, min_lat, max_lat
    real(r_kind) x1(3), x2(3), x3(3), x4(3), x(3)
    integer is1, ie1, is2, ie2, js, je
    real(r_kind) :: q(-rw:rw,-rw:rw,1)
    real(r_kind) :: fo(neq*4+ncq)
    real(r_kind) :: slmsk_threshold = 0

    call log_notice('Calculating averaged surface height and land fraction on model grids.', pid=proc%id)
    topo_dlon = topo_lon(2) - topo_lon(1)
    topo_dlat = topo_lat(2) - topo_lat(1)
    
    associate (mesh => static%mesh)
      do idom = 1, mesh%ndom
        if (.not. mesh%domain(idom)%initialized) cycle
        associate (domain   => mesh%domain(idom),            &
                   lon      => mesh%domain(idom)%lon,        &
                   lat      => mesh%domain(idom)%lat,        &
                   zs       => static%domain(idom)%zs,       &
                   landfrac => static%domain(idom)%landfrac, &
                   slmsk    => static%domain(idom)%slmsk)
        call mesh%get_params(idom=idom, pvs=pvs, pve=pve, irs=irs, ire=ire, jrs=jrs, jre=jre, ims=ims, ime=ime, jms=jms, jme=jme)
        !$OMP PARALLEL DO PRIVATE(min_lon,max_lon,min_lat,max_lat,x1,x2,x3,x4,x,is1,ie1,is2,ie2,js,je,n,i,j) COLLAPSE(2)
        do jg = jms, jme
          do ig = ims, ime
            min_lon = minval(lon(pvs:pve,ig,jg))
            max_lon = maxval(lon(pvs:pve,ig,jg))
            min_lat = minval(lat(pvs:pve,ig,jg))
            max_lat = maxval(lat(pvs:pve,ig,jg))
      
            call latlon_to_cart(lon(pvs+0,ig,jg), lat(2,ig,jg), x1(1), x1(2), x1(3))
            call latlon_to_cart(lon(pvs+1,ig,jg), lat(3,ig,jg), x2(1), x2(2), x2(3))
            call latlon_to_cart(lon(pvs+2,ig,jg), lat(4,ig,jg), x3(1), x3(2), x3(3))
            call latlon_to_cart(lon(pvs+3,ig,jg), lat(5,ig,jg), x4(1), x4(2), x4(3))
      
            if (max_lon - min_lon > 300 * D2R) then
              is1 = int(max_lon / topo_dlon) + 1
              ie1 = topo_nlon
              is2 = 1
              ie2 = int(min_lon / topo_dlon) + 1
              if (ie2 == 1) ie2 = 0 ! ETOPO1 has 0 and 360 longitude, so drop one.
            else
              is1 = int(min_lon / topo_dlon) + 1
              ie1 = min( int(max_lon / topo_dlon) + 1, topo_nlon )
              is2 = 1
              ie2 = merge( 1, 0, ie1==topo_nlon )
            end if
            js = int((min_lat + pi0p5) / topo_dlat) + 1
            je = int((max_lat + pi0p5) / topo_dlat) + 1
      
            zs(ig,jg) = 0
            if(allocated(topo_landfrac))landfrac(ig,jg) = 0
            n = 0
            do j = js, je
              do i = is1, ie1
                call latlon_to_cart(topo_lon(i), topo_lat(j), x(1), x(2), x(3))
                if (orient(x1, x2, x) == 1) cycle
                if (orient(x2, x3, x) == 1) cycle
                if (orient(x3, x4, x) == 1) cycle
                if (orient(x4, x1, x) == 1) cycle
                zs(ig,jg) = zs(ig,jg) + topo_zs(i,j)
                if(allocated(topo_landfrac))then
                  if(topo_landfrac(i,j)>slmsk_threshold)landfrac(ig,jg) = landfrac(ig,jg) + 1
                endif
                
                n = n + 1
              end do
              do i = is2, ie2
                call latlon_to_cart(topo_lon(i), topo_lat(j), x(1), x(2), x(3))
                if (orient(x1, x2, x) == 1) cycle
                if (orient(x2, x3, x) == 1) cycle
                if (orient(x3, x4, x) == 1) cycle
                if (orient(x4, x1, x) == 1) cycle
                zs(ig,jg) = zs(ig,jg) + topo_zs(i,j)
                if(allocated(topo_landfrac))then
                  if(topo_landfrac(i,j)>slmsk_threshold)landfrac(ig,jg) = landfrac(ig,jg) + 1
                endif
                
                n = n + 1
              end do
            end do
            
            zs(ig,jg) = zs(ig,jg) / n
            if(allocated(topo_landfrac))then
              landfrac(ig,jg) = landfrac(ig,jg) / n
              slmsk(ig,jg) = 1
              if( landfrac(ig,jg)<=slmsk_threshold .or. zs(ig,jg)<0 )then
                zs   (ig,jg) = 0 ! Remove water grid
                slmsk(ig,jg) = 0
              endif
            endif
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end do
    end associate
    
    call fill_halo(static%array)
    
    call topo_filter(static)
    
  end subroutine topo_regrid
  
  subroutine topo_filter(static)
    type(static_type), intent(inout) :: static

    integer idom, ip, p, i, j, ig, jg, n, irow, icol
    integer pc, pqs, pqe, pes(6), pee(6), ims, ime, jms, jme, its, ite, jts, jte
    integer irs, ire, jrs, jre
    real(r_kind) min_lon, max_lon, min_lat, max_lat
    real(r_kind) x1(3), x2(3), x3(3), x4(3), x(3)
    integer is1, ie1, is2, ie2, js, je
    real(r_kind) :: q(-rw:rw,-rw:rw,1)
    real(r_kind) :: fo(neq*4+ncq)
    real(r_kind), dimension(:,:), allocatable :: zs_tmp

    real(r_kind) :: JhC, JhCL, JhCR, JhCB, JhCT
    real(r_kind) :: JhL, JhR, JhB, JhT
    real(r_kind), dimension(2,2) :: iGhL, iGhR, iGhB, iGhT
    real(r_kind), dimension(2,2) :: iGhCL, iGhCR, iGhCB, iGhCT
    real(r_kind) lap, lap_xL, lap_xR, lap_xB, lap_xT, lap_yL, lap_yR, lap_yB, lap_yT, lap_xx, lap_xy, lap_yy, lap_yx
    
    integer :: loop
    integer :: filt_loop = 0
    
    do loop = 1, filt_loop
      associate (mesh => static%mesh)
        do idom = 1, mesh%ndom
          if (.not. mesh%domain(idom)%initialized) cycle
          associate (domain   => mesh%domain(idom),            &
                     dx       => mesh%domain(idom)%dx,         &
                     dy       => mesh%domain(idom)%dy,         &
                     lon      => mesh%domain(idom)%lon,        &
                     lat      => mesh%domain(idom)%lat,        &
                     Jh       => mesh%domain(idom)%Jh,         &
                     iGh      => mesh%domain(idom)%iGh,        &
                     zs       => static%domain(idom)%zs,       &
                     landfrac => static%domain(idom)%landfrac, &
                     slmsk    => static%domain(idom)%slmsk     )
          call mesh%get_params(idom=idom, its=its, ite=ite, jts=jts, jte=jte, ims=ims, ime=ime, jms=jms, jme=jme, pc=pc, pqs=pqs, pqe=pqe, pes=pes, pee=pee)

          if(allocated(zs_tmp))deallocate(zs_tmp)
          allocate(zs_tmp(ims:ime,jms:jme))
          zs_tmp = zs
          
          !!$OMP PARALLEL DO COLLAPSE(2)
          !do j = jts, jte
          !  do i = its, ite
          !    !if( slmsk(i,j)==1 )then
          !      zs(i,j) = zs_tmp(i,j) + 0.125 * ( zs_tmp(i+1,j) + zs_tmp(i-1,j) + zs_tmp(i,j+1) + zs_tmp(i,j-1) - 4 * zs_tmp(i,j) )
          !      !zs(i,j) = zs_tmp(i,j) - 0.03125 * ( zs_tmp(i+2,j) + zs_tmp(i-2,j) + zs_tmp(i,j+2) + zs_tmp(i,j-2) - 4 * ( zs_tmp(i+1,j) + zs_tmp(i-1,j) + zs_tmp(i,j+1) + zs_tmp(i,j-1)) + 12 * zs_tmp(i,j) )
          !      !zs(i,j) = zs_tmp(i,j) + 0.0078125*(zs_tmp(i+3,j)+zs_tmp(i-3,j)+zs_tmp(i,j+3)+zs_tmp(i,j-3)-6*(zs_tmp(i+2,j)+zs_tmp(i-2,j)+zs_tmp(i,j+2)+zs_tmp(i,j-2))+15*(zs_tmp(i+1,j)+zs_tmp(i-1,j)+zs_tmp(i,j+1)+zs_tmp(i,j-1))-40*zs_tmp(i,j))
          !      !zs(i,j) = zs_tmp(i,j) - 0.001953125*(zs_tmp(i+4,j)+zs_tmp(i-4,j)+zs_tmp(i,j+4)+zs_tmp(i,j-4)-8*(zs_tmp(i+3,j)+zs_tmp(i-3,j)+zs_tmp(i,j+3)+zs_tmp(i,j-3))+28*(zs_tmp(i+2,j)+zs_tmp(i-2,j)+zs_tmp(i,j+2)+zs_tmp(i,j-2))-56*(zs_tmp(i+1,j)+zs_tmp(i-1,j)+zs_tmp(i,j+1)+zs_tmp(i,j-1))+140*zs_tmp(i,j))
          !    !endif
          !  end do
          !end do
          !!$OMP END PARALLEL DO
          
          !$OMP PARALLEL DO PRIVATE(ip, p, JhC, JhCL, JhCR, JhCB, JhCT, JhL, JhR, JhB, JhT, &
          !$OMP                     iGhL, iGhR, iGhB, iGhT, iGhCL, iGhCR, iGhCB, iGhCT, &
          !$OMP                     lap, lap_xL, lap_xR, lap_xB, lap_xT, lap_yL, lap_yR, lap_yB, lap_yT, &
          !$OMP                     lap_xx, lap_xy, lap_yy, lap_yx) COLLAPSE(2)
          do j = jts, jte
            do i = its, ite
              !if( slmsk(i,j)==1 )then
                JhC  = mesh%domain(idom)%cell_quad(Jh(pqs:pqe,i  ,j  ))
                JhCL = mesh%domain(idom)%cell_quad(Jh(pqs:pqe,i-1,j  ))
                JhCR = mesh%domain(idom)%cell_quad(Jh(pqs:pqe,i+1,j  ))
                JhCB = mesh%domain(idom)%cell_quad(Jh(pqs:pqe,i  ,j-1))
                JhCT = mesh%domain(idom)%cell_quad(Jh(pqs:pqe,i  ,j+1))
                
                JhL = mesh%domain(idom)%edge_quad(Jh(pes(1):pee(1),i,j))
                JhR = mesh%domain(idom)%edge_quad(Jh(pes(2):pee(2),i,j))
                JhB = mesh%domain(idom)%edge_quad(Jh(pes(3):pee(3),i,j))
                JhT = mesh%domain(idom)%edge_quad(Jh(pes(4):pee(4),i,j))
                
                do icol = 1, 2
                  do irow = 1, 2
                    iGhCL(irow,icol) = mesh%domain(idom)%cell_quad(iGh(irow,icol,pqs:pqe,i-1,j  ))
                    iGhCR(irow,icol) = mesh%domain(idom)%cell_quad(iGh(irow,icol,pqs:pqe,i+1,j  ))
                    iGhCB(irow,icol) = mesh%domain(idom)%cell_quad(iGh(irow,icol,pqs:pqe,i  ,j-1))
                    iGhCT(irow,icol) = mesh%domain(idom)%cell_quad(iGh(irow,icol,pqs:pqe,i  ,j+1))
                    
                    iGhL(irow,icol) = mesh%domain(idom)%edge_quad(iGh(irow,icol,pes(1):pee(1),i,j))
                    iGhR(irow,icol) = mesh%domain(idom)%edge_quad(iGh(irow,icol,pes(2):pee(2),i,j))
                    iGhB(irow,icol) = mesh%domain(idom)%edge_quad(iGh(irow,icol,pes(3):pee(3),i,j))
                    iGhT(irow,icol) = mesh%domain(idom)%edge_quad(iGh(irow,icol,pes(4):pee(4),i,j))
                  enddo
                enddo
                
                lap_xL = JhR  * iGhR (1,1) * (zs_tmp(i  ,j  ) - zs_tmp(i-1,j  )) / dx 
                lap_xR = JhL  * iGhL (1,1) * (zs_tmp(i+1,j  ) - zs_tmp(i  ,j  )) / dx 
                lap_xB = JhCB * iGhCB(2,1) * (zs_tmp(i+1,j-1) - zs_tmp(i-1,j-1)) / ( 2 * dy ) 
                lap_xT = JhCT * iGhCT(2,1) * (zs_tmp(i+1,j+1) - zs_tmp(i-1,j+1)) / ( 2 * dy ) 
                lap_yL = JhCL * iGhCR(1,2) * (zs_tmp(i-1,j+1) - zs_tmp(i-1,j-1)) / ( 2 * dx ) 
                lap_yR = JhCR * iGhCL(1,2) * (zs_tmp(i+1,j+1) - zs_tmp(i+1,j-1)) / ( 2 * dx ) 
                lap_yB = JhB  * iGhB (2,2) * (zs_tmp(i  ,j  ) - zs_tmp(i  ,j-1)) / dy 
                lap_yT = JhT  * iGhT (2,2) * (zs_tmp(i  ,j+1) - zs_tmp(i  ,j  )) / dy 
                
                
                lap_xx = ( lap_xR - lap_xL ) / dx
                lap_xy = ( lap_yR - lap_yL ) / ( 2 * dx )
                lap_yy = ( lap_yT - lap_yB ) / dy
                lap_yx = ( lap_xT - lap_xT ) / ( 2 * dy )
                
                !lap_xx = ( lap_xR * min( 1._r_kind, slmsk(i,j), slmsk(i+1,j) ) - lap_xL * min( 1._r_kind, slmsk(i-1,j), slmsk(i,j) ) ) / dx
                !lap_xy = ( lap_yR * min( 1._r_kind, slmsk(i,j), slmsk(i+1,j) ) - lap_yL * min( 1._r_kind, slmsk(i-1,j), slmsk(i,j) ) ) / ( 2 * dx )
                !lap_yy = ( lap_yT * min( 1._r_kind, slmsk(i,j), slmsk(i,j+1) ) - lap_yB * min( 1._r_kind, slmsk(i,j-1), slmsk(i,j) ) ) / dy
                !lap_yx = ( lap_xT * min( 1._r_kind, slmsk(i,j), slmsk(i,j+1) ) - lap_xT * min( 1._r_kind, slmsk(i,j-1), slmsk(i,j) ) ) / ( 2 * dy )
                
                lap = ( lap_xx + lap_xy + lap_yy + lap_yx ) / JhC
                
                zs(i,j) = zs_tmp(i,j) + 0.125 * dx * dy * lap
              !endif
            end do
          end do
          !$OMP END PARALLEL DO
          end associate
        enddo
      end associate
      call fill_halo(static%array)
      call ghost_interp_topo(static)
    enddo
    
  end subroutine topo_filter
  
  subroutine topo_recon(static)
    type(static_type), intent(inout) :: static

    integer ims, ime, jms, jme, kms, kme
    integer ids, ide, jds, jde, kds, kde
    integer irs, ire, jrs, jre
    integer igs, ige, jgs, jge
    integer neq, nfq, nfqv
    integer idom, p, i, j, k
    integer idom_ngb, ps, is, js
    integer pqs, pqe
  
    call log_notice('Calculating topo reconstruction.', pid=proc%id)
    
    call reconstruct_topo(static)
    
    associate (mesh => static%mesh)
      do idom = 1, mesh%ndom
        if (.not. mesh%domain(idom)%initialized) cycle
        call mesh%get_params(idom=idom, neq=neq, nfq=nfq, nfqv=nfqv, &
                             ims=ims, ime=ime, jms=jms, jme=jme,     &
                             ids=ids, ide=ide, jds=jds, jde=jde,     &
                             kms=kms, kme=kme, kds=kds, kde=kde,     &
                             irs=irs, ire=ire, jrs=jrs, jre=jre,     &
                             igs=igs, ige=ige, jgs=jgs, jge=jge,     &
                             pqs=pqs, pqe=pqe                        )
        associate (dx     => mesh%domain(idom)%dx,      &
                   dy     => mesh%domain(idom)%dy,      &
                   slmsk  => static%domain(idom)%slmsk, &
                   zs     => static%domain(idom)%zs,    &
                   zsC    => static%domain(idom)%zsC,   &
                   zsL    => static%domain(idom)%zsL,   &
                   zsR    => static%domain(idom)%zsR,   &
                   zsB    => static%domain(idom)%zsB,   &
                   zsT    => static%domain(idom)%zsT,   &
                   zsD    => static%domain(idom)%zsD,   &
                   zsU    => static%domain(idom)%zsU,   &
                   zsQ    => static%domain(idom)%zsQ    )
        !$OMP PARALLEL DO COLLAPSE(2)
        do j = jms, jme
          do i = ims, ime
            zsC(i,j) = zs(i,j)
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end do
    end associate
  end subroutine topo_recon
  
  subroutine topo_recon_deriv(static)
    type(static_type), intent(inout) :: static

    integer ims, ime, jms, jme, kms, kme
    integer ids, ide, jds, jde, kds, kde
    integer irs, ire, jrs, jre
    integer its, ite, jts, jte
    integer igs, ige, jgs, jge
    integer neq, nfq, nfqv
    integer idom, p, i, j, k
    integer idom_ngb, ps, is, js
    integer pqs, pqe
    integer pes(6), pee(6)
    real(r_kind) :: msk
  
    call log_notice('Calculating topo derivative.', pid=proc%id)
    
    call reconstruct_topo_deriv(static)
    
    associate (mesh => static%mesh)
      do idom = 1, mesh%ndom
        if (.not. mesh%domain(idom)%initialized) cycle
        call mesh%get_params(idom=idom, neq=neq, nfq=nfq, nfqv=nfqv, &
                             ims=ims, ime=ime, jms=jms, jme=jme,     &
                             ids=ids, ide=ide, jds=jds, jde=jde,     &
                             kms=kms, kme=kme, kds=kds, kde=kde,     &
                             irs=irs, ire=ire, jrs=jrs, jre=jre,     &
                             its=its, ite=ite, jts=jts, jte=jte,     &
                             igs=igs, ige=ige, jgs=jgs, jge=jge,     &
                             pqs=pqs, pqe=pqe, pes=pes, pee=pee      )
        associate (dx     => mesh%domain(idom)%dx,      &
                   dy     => mesh%domain(idom)%dy,      &
                   zs     => static%domain(idom)%zs,    &
                   zsC    => static%domain(idom)%zsC,   &
                   zsL    => static%domain(idom)%zsL,   &
                   zsR    => static%domain(idom)%zsR,   &
                   zsB    => static%domain(idom)%zsB,   &
                   zsT    => static%domain(idom)%zsT,   &
                   zsD    => static%domain(idom)%zsD,   &
                   zsU    => static%domain(idom)%zsU,   &
                   zsQ    => static%domain(idom)%zsQ,   &    
                   slmsk  => static%domain(idom)%slmsk, &
                   dzsdx  => static%domain(idom)%dzsdx, &
                   dzsdxC => static%domain(idom)%dzsdxC,&
                   dzsdxL => static%domain(idom)%dzsdxL,&
                   dzsdxR => static%domain(idom)%dzsdxR,&
                   dzsdxB => static%domain(idom)%dzsdxB,&
                   dzsdxT => static%domain(idom)%dzsdxT,&
                   dzsdxD => static%domain(idom)%dzsdxD,&
                   dzsdxU => static%domain(idom)%dzsdxU,&
                   dzsdxQ => static%domain(idom)%dzsdxQ,&
                   dzsdy  => static%domain(idom)%dzsdy, &
                   dzsdyC => static%domain(idom)%dzsdyC,&
                   dzsdyL => static%domain(idom)%dzsdyL,&
                   dzsdyR => static%domain(idom)%dzsdyR,&
                   dzsdyB => static%domain(idom)%dzsdyB,&
                   dzsdyT => static%domain(idom)%dzsdyT,&
                   dzsdyD => static%domain(idom)%dzsdyD,&
                   dzsdyU => static%domain(idom)%dzsdyU,&
                   dzsdyQ => static%domain(idom)%dzsdyQ )
        !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(p)
        do j = jms, jme
          do i = ims, ime
            
            dzsdx(i,j) = mesh%domain(idom)%cell_quad( dzsdxQ(:,i,j) )
            dzsdy(i,j) = mesh%domain(idom)%cell_quad( dzsdyQ(:,i,j) )
            
            dzsdxC(i,j) = dzsdx(i,j)
            dzsdyC(i,j) = dzsdy(i,j)
            
            !if( i>ims+3 .and. i<ime-3 .and. j>jms+3 .and. j<jme-3 )then
            !  msk = 1!min(slmsk(i,j),1.) !min( sum(slmsk(i-1:i+1,j-1:j+1)), 1. )
            !
            !  !dzsdx(i,j) = center_diff2(zs(i-1:i+1,j), dx)
            !  !dzsdy(i,j) = center_diff2(zs(i,j-1:j+1), dy)
            !  !
            !  !dzsdxC(i,j) = center_diff2(zsC(i-1:i+1,j), dx)
            !  !dzsdyC(i,j) = center_diff2(zsC(i,j-1:j+1), dy)
            !  !
            !  !do p = pes(1),pee(1)
            !  !  dzsdxL(p,i,j) = center_diff2(zsL(p,i-1:i+1,j), dx)
            !  !  dzsdyL(p,i,j) = center_diff2(zsL(p,i,j-1:j+1), dy)
            !  !enddo
            !  !do p = pes(2),pee(2)
            !  !  dzsdxR(p,i,j) = center_diff2(zsR(p,i-1:i+1,j), dx)
            !  !  dzsdyR(p,i,j) = center_diff2(zsR(p,i,j-1:j+1), dy)
            !  !enddo
            !  !do p = pes(3),pee(3)
            !  !  dzsdxB(p,i,j) = center_diff2(zsB(p,i-1:i+1,j), dx)
            !  !  dzsdyB(p,i,j) = center_diff2(zsB(p,i,j-1:j+1), dy)
            !  !enddo
            !  !do p = pes(4),pee(4)
            !  !  dzsdxT(p,i,j) = center_diff2(zsT(p,i-1:i+1,j), dx)
            !  !  dzsdyT(p,i,j) = center_diff2(zsT(p,i,j-1:j+1), dy)
            !  !enddo
            !  !do p = pqs,pqe
            !  !  dzsdxQ(p,i,j) = center_diff2(zsQ(p,i-1:i+1,j), dx)
            !  !  dzsdyQ(p,i,j) = center_diff2(zsQ(p,i,j-1:j+1), dy)
            !  !enddo
            !  
            !  !dzsdx(i,j) = center_diff4(zs(i-2:i+2,j), dx)
            !  !dzsdy(i,j) = center_diff4(zs(i,j-2:j+2), dy)
            !  !
            !  !dzsdxC(i,j) = center_diff4(zsC(i-2:i+2,j), dx)
            !  !dzsdyC(i,j) = center_diff4(zsC(i,j-2:j+2), dy)
            !  !
            !  !do p = pes(1),pee(1)
            !  !  dzsdxL(p,i,j) = center_diff4(zsL(p,i-2:i+2,j), dx)
            !  !  dzsdyL(p,i,j) = center_diff4(zsL(p,i,j-2:j+2), dy)
            !  !enddo
            !  !do p = pes(2),pee(2)
            !  !  dzsdxR(p,i,j) = center_diff4(zsR(p,i-2:i+2,j), dx)
            !  !  dzsdyR(p,i,j) = center_diff4(zsR(p,i,j-2:j+2), dy)
            !  !enddo
            !  !do p = pes(3),pee(3)
            !  !  dzsdxB(p,i,j) = center_diff4(zsB(p,i-2:i+2,j), dx)
            !  !  dzsdyB(p,i,j) = center_diff4(zsB(p,i,j-2:j+2), dy)
            !  !enddo
            !  !do p = pes(4),pee(4)
            !  !  dzsdxT(p,i,j) = center_diff4(zsT(p,i-2:i+2,j), dx)
            !  !  dzsdyT(p,i,j) = center_diff4(zsT(p,i,j-2:j+2), dy)
            !  !enddo
            !  !do p = pqs,pqe
            !  !  dzsdxQ(p,i,j) = center_diff4(zsQ(p,i-2:i+2,j), dx)
            !  !  dzsdyQ(p,i,j) = center_diff4(zsQ(p,i,j-2:j+2), dy)
            !  !enddo
            !  
            !  !dzsdx(i,j) = center_diff6(zs(i-3:i+3,j), dx)
            !  !dzsdy(i,j) = center_diff6(zs(i,j-3:j+3), dy)
            !  !
            !  !dzsdxC(i,j) = center_diff6(zsC(i-3:i+3,j), dx)
            !  !dzsdyC(i,j) = center_diff6(zsC(i,j-3:j+3), dy)
            !  !
            !  !do p = pes(1),pee(1)
            !  !  dzsdxL(p,i,j) = center_diff6(zsL(p,i-3:i+3,j), dx)
            !  !  dzsdyL(p,i,j) = center_diff6(zsL(p,i,j-3:j+3), dy)
            !  !enddo
            !  !do p = pes(2),pee(2)
            !  !  dzsdxR(p,i,j) = center_diff6(zsR(p,i-3:i+3,j), dx)
            !  !  dzsdyR(p,i,j) = center_diff6(zsR(p,i,j-3:j+3), dy)
            !  !enddo
            !  !do p = pes(3),pee(3)
            !  !  dzsdxB(p,i,j) = center_diff6(zsB(p,i-3:i+3,j), dx)
            !  !  dzsdyB(p,i,j) = center_diff6(zsB(p,i,j-3:j+3), dy)
            !  !enddo
            !  !do p = pes(4),pee(4)
            !  !  dzsdxT(p,i,j) = center_diff6(zsT(p,i-3:i+3,j), dx)
            !  !  dzsdyT(p,i,j) = center_diff6(zsT(p,i,j-3:j+3), dy)
            !  !enddo
            !  !do p = pqs,pqe
            !  !  dzsdxQ(p,i,j) = center_diff6(zsQ(p,i-3:i+3,j), dx)
            !  !  dzsdyQ(p,i,j) = center_diff6(zsQ(p,i,j-3:j+3), dy)
            !  !enddo
            !
            !  ! 8th order
            !  dzsdx(i,j) = center_diff8(zs(i-4:i+4,j), dx)
            !  dzsdy(i,j) = center_diff8(zs(i,j-4:j+4), dy)
            !  
            !  dzsdxC(i,j) = center_diff8(zsC(i-4:i+4,j), dx)
            !  dzsdyC(i,j) = center_diff8(zsC(i,j-4:j+4), dy)
            !  
            !  do p = pes(1),pee(1)
            !    dzsdxL(p,i,j) = center_diff8(zsL(p,i-4:i+4,j), dx)
            !    dzsdyL(p,i,j) = center_diff8(zsL(p,i,j-4:j+4), dy)
            !  enddo
            !  do p = pes(2),pee(2)
            !    dzsdxR(p,i,j) = center_diff8(zsR(p,i-4:i+4,j), dx)
            !    dzsdyR(p,i,j) = center_diff8(zsR(p,i,j-4:j+4), dy)
            !  enddo
            !  do p = pes(3),pee(3)
            !    dzsdxB(p,i,j) = center_diff8(zsB(p,i-4:i+4,j), dx)
            !    dzsdyB(p,i,j) = center_diff8(zsB(p,i,j-4:j+4), dy)
            !  enddo
            !  do p = pes(4),pee(4)
            !    dzsdxT(p,i,j) = center_diff8(zsT(p,i-4:i+4,j), dx)
            !    dzsdyT(p,i,j) = center_diff8(zsT(p,i,j-4:j+4), dy)
            !  enddo
            !  do p = pqs,pqe
            !    dzsdxQ(p,i,j) = center_diff8(zsQ(p,i-4:i+4,j), dx)
            !    dzsdyQ(p,i,j) = center_diff8(zsQ(p,i,j-4:j+4), dy)
            !  enddo
            !
            !  !! 12th order
            !  !dzsdx(i,j) = center_diff12(zs(i-6:i+6,j), dx)
            !  !dzsdy(i,j) = center_diff12(zs(i,j-6:j+6), dy)
            !  !
            !  !dzsdxC(i,j) = center_diff12(zsC(i-6:i+6,j), dx)
            !  !dzsdyC(i,j) = center_diff12(zsC(i,j-6:j+6), dy)
            !  !
            !  !do p = pes(1),pee(1)
            !  !  dzsdxL(p,i,j) = center_diff12(zsL(p,i-6:i+6,j), dx)
            !  !  dzsdyL(p,i,j) = center_diff12(zsL(p,i,j-6:j+6), dy)
            !  !enddo
            !  !do p = pes(2),pee(2)
            !  !  dzsdxR(p,i,j) = center_diff12(zsR(p,i-6:i+6,j), dx)
            !  !  dzsdyR(p,i,j) = center_diff12(zsR(p,i,j-6:j+6), dy)
            !  !enddo
            !  !do p = pes(3),pee(3)
            !  !  dzsdxB(p,i,j) = center_diff12(zsB(p,i-6:i+6,j), dx)
            !  !  dzsdyB(p,i,j) = center_diff12(zsB(p,i,j-6:j+6), dy)
            !  !enddo
            !  !do p = pes(4),pee(4)
            !  !  dzsdxT(p,i,j) = center_diff12(zsT(p,i-6:i+6,j), dx)
            !  !  dzsdyT(p,i,j) = center_diff12(zsT(p,i,j-6:j+6), dy)
            !  !enddo
            !  !do p = pqs,pqe
            !  !  dzsdxQ(p,i,j) = center_diff12(zsQ(p,i-6:i+6,j), dx)
            !  !  dzsdyQ(p,i,j) = center_diff12(zsQ(p,i,j-6:j+6), dy)
            !  !enddo
            !
            !  dzsdx(i,j) = dzsdx(i,j) * msk
            !  dzsdy(i,j) = dzsdy(i,j) * msk
            !  
            !  dzsdxC(i,j) = dzsdxC(i,j) * msk
            !  dzsdyC(i,j) = dzsdyC(i,j) * msk
            !  
            !  do p = pes(1),pee(1)
            !    dzsdxL(p,i,j) = dzsdxL(p,i,j) * msk
            !    dzsdyL(p,i,j) = dzsdyL(p,i,j) * msk
            !  enddo
            !  do p = pes(2),pee(2)
            !    dzsdxR(p,i,j) = dzsdxR(p,i,j) * msk
            !    dzsdyR(p,i,j) = dzsdyR(p,i,j) * msk
            !  enddo
            !  do p = pes(3),pee(3)
            !    dzsdxB(p,i,j) = dzsdxB(p,i,j) * msk
            !    dzsdyB(p,i,j) = dzsdyB(p,i,j) * msk
            !  enddo
            !  do p = pes(4),pee(4)
            !    dzsdxT(p,i,j) = dzsdxT(p,i,j) * msk
            !    dzsdyT(p,i,j) = dzsdyT(p,i,j) * msk
            !  enddo
            !  do p = pqs,pqe
            !    dzsdxQ(p,i,j) = dzsdxQ(p,i,j) * msk
            !    dzsdyQ(p,i,j) = dzsdyQ(p,i,j) * msk
            !  enddo
            !
            !
            !endif
          end do
        end do
        !$OMP END PARALLEL DO
        end associate
      end do
    end associate
  end subroutine topo_recon_deriv

  subroutine topo_unify_edges(static)
    type(static_type), intent(inout) :: static

    integer idom, idom_ngb, orient, orient_adj
    integer i, j, k, p, i_ngb, j_ngb, p_ngb, i_adj, j_adj, p_adj
    integer ids, ide, jds, jde, kds, kde, pes(6), pee(6)
    integer is, ie, js, je, di, dj
    real(r8), pointer, dimension(:,:,:) :: zs, zs_ngb
    
    call log_notice('Unifying topo on cell edges.', pid=proc%id)

    associate (mesh => static%mesh)
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%active) then
          call mesh%get_params(idom=idom, ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, pes=pes, pee=pee)
          do orient = 1, 4
            if (mesh%domain(idom)%touch_edges(orient)) then
              idom_ngb = mesh%domain(idom)%ngb(orient)%ptr%id
              associate (ngb => mesh%domain(idom)%ngb(orient))
              is = ids; ie = ide; js = jds; je = jde; di = 0; dj = 0
              if (orient == 1) then; ie = ids; di = -1; orient_adj = 2; zs => static%domain(idom)%zsR; end if
              if (orient == 2) then; is = ide; di = +1; orient_adj = 1; zs => static%domain(idom)%zsL; end if
              if (orient == 3) then; je = jds; dj = -1; orient_adj = 4; zs => static%domain(idom)%zsT; end if
              if (orient == 4) then; js = jde; dj = +1; orient_adj = 3; zs => static%domain(idom)%zsB; end if
              if (ngb%orient == 1) zs_ngb => static%domain(idom_ngb)%zsL
              if (ngb%orient == 2) zs_ngb => static%domain(idom_ngb)%zsR
              if (ngb%orient == 3) zs_ngb => static%domain(idom_ngb)%zsB
              if (ngb%orient == 4) zs_ngb => static%domain(idom_ngb)%zsT
              !$OMP PARALLEL DO PRIVATE(p, i_ngb, j_ngb, p_ngb, i_adj, j_adj, p_adj) COLLAPSE(3)
              do j = js, je
                do i = is, ie
                  do p = pes(orient), pee(orient)
                    i_adj = i + di; j_adj = j + dj
                    p_adj = p - pes(orient) + pes(orient_adj)
                    call mesh%translate_idx_along_edge(idom, orient, i, j, p, i_ngb, j_ngb, p_ngb)
                    
                    zs(p_adj,i_adj,j_adj) = zs_ngb(p_ngb,i_ngb,j_ngb)
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            end if
          end do
        end if
      end do
      
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%active) then
          call mesh%get_params(idom=idom, ids=ids, ide=ide, jds=jds, jde=jde)
          !$OMP PARALLEL
          !$OMP DO COLLAPSE(2)
          do j = jds, jde
            do i = ids, ide + 1
              ! unify zs on x-dir
              static%domain(idom)%zsL(:,i,j) = 0.5 * ( static%domain(idom)%zsL(:,i,j) + static%domain(idom)%zsR(:,i-1,j) )
              static%domain(idom)%zsR(:,i-1,j) = static%domain(idom)%zsL(:,i,j)
            end do
          end do
          !$OMP END DO
          !$OMP DO COLLAPSE(2)
          do j = jds, jde + 1
            do i = ids, ide
              ! unify zs on y-dir
              static%domain(idom)%zsB(:,i,j) = 0.5 * ( static%domain(idom)%zsB(:,i,j) + static%domain(idom)%zsT(:,i,j-1) )
              static%domain(idom)%zsT(:,i,j-1) = static%domain(idom)%zsB(:,i,j)
            end do
          end do
          !$OMP END DO
          !$OMP END PARALLEL
        end if
      end do
    end associate
    
  end subroutine topo_unify_edges

  subroutine topo_unify_edges_deriv(static)
    type(static_type), intent(inout) :: static

    integer idom, idom_ngb, orient, orient_adj
    integer i, j, k, p, i_ngb, j_ngb, p_ngb, i_adj, j_adj, p_adj
    integer ids, ide, jds, jde, kds, kde, pes(6), pee(6)
    integer is, ie, js, je, di, dj
    real(r8), pointer, dimension(:,:,:) :: dzsdx, dzsdx_ngb
    real(r8), pointer, dimension(:,:,:) :: dzsdy, dzsdy_ngb
    
    real(r_kind) dzsdlon, dzsdlat
    real(r_kind) Ah(2,2), iAh(2,2)
    
    call log_notice('Unifying topo derivatve on cell edges.', pid=proc%id)

    associate (mesh => static%mesh)
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%active) then
          call mesh%get_params(idom=idom, ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, pes=pes, pee=pee)
          do orient = 1, 4
            if (mesh%domain(idom)%touch_edges(orient)) then
              idom_ngb = mesh%domain(idom)%ngb(orient)%ptr%id
              associate (ngb    => mesh%domain(idom)%ngb(orient))
              is = ids; ie = ide; js = jds; je = jde; di = 0; dj = 0
              if (orient == 1) then; ie = ids; di = -1; orient_adj = 2; dzsdx => static%domain(idom)%dzsdxR; end if
              if (orient == 2) then; is = ide; di = +1; orient_adj = 1; dzsdx => static%domain(idom)%dzsdxL; end if
              if (orient == 3) then; je = jds; dj = -1; orient_adj = 4; dzsdx => static%domain(idom)%dzsdxT; end if
              if (orient == 4) then; js = jde; dj = +1; orient_adj = 3; dzsdx => static%domain(idom)%dzsdxB; end if
              if (orient == 1) then; ie = ids; di = -1; orient_adj = 2; dzsdy => static%domain(idom)%dzsdyR; end if
              if (orient == 2) then; is = ide; di = +1; orient_adj = 1; dzsdy => static%domain(idom)%dzsdyL; end if
              if (orient == 3) then; je = jds; dj = -1; orient_adj = 4; dzsdy => static%domain(idom)%dzsdyT; end if
              if (orient == 4) then; js = jde; dj = +1; orient_adj = 3; dzsdy => static%domain(idom)%dzsdyB; end if
              if (ngb%orient == 1) dzsdx_ngb => static%domain(idom_ngb)%dzsdxL
              if (ngb%orient == 2) dzsdx_ngb => static%domain(idom_ngb)%dzsdxR
              if (ngb%orient == 3) dzsdx_ngb => static%domain(idom_ngb)%dzsdxB
              if (ngb%orient == 4) dzsdx_ngb => static%domain(idom_ngb)%dzsdxT
              if (ngb%orient == 1) dzsdy_ngb => static%domain(idom_ngb)%dzsdyL
              if (ngb%orient == 2) dzsdy_ngb => static%domain(idom_ngb)%dzsdyR
              if (ngb%orient == 3) dzsdy_ngb => static%domain(idom_ngb)%dzsdyB
              if (ngb%orient == 4) dzsdy_ngb => static%domain(idom_ngb)%dzsdyT
              k = 1
              !$OMP PARALLEL DO PRIVATE(p, i_ngb, j_ngb, p_ngb, i_adj, j_adj, p_adj, Ah, iAh) COLLAPSE(3)
              do j = js, je
                do i = is, ie
                  do p = pes(orient), pee(orient)
                    i_adj = i + di; j_adj = j + dj
                    p_adj = p - pes(orient) + pes(orient_adj)
                    call mesh%translate_idx_along_edge(idom, orient, i, j, p, i_ngb, j_ngb, p_ngb)
                    
                    iAh = mesh%get_dom_iAh(idom_ngb,p_ngb,i_ngb,j_ngb,k)
                    call convert_hor_deriv_cube_to_sph(dzsdlon,dzsdlat,dzsdx_ngb(p_ngb,i_ngb,j_ngb),dzsdy_ngb(p_ngb,i_ngb,j_ngb),iAh)
                    
                    Ah = mesh%get_dom_Ah(idom,p_adj,i_adj,j_adj,k)
                    call convert_hor_deriv_sph_to_cube(dzsdx(p_adj,i_adj,j_adj),dzsdy(p_adj,i_adj,j_adj),dzsdlon,dzsdlat,Ah)
                  end do
                end do
              end do
              !$OMP END PARALLEL DO
              end associate
            end if
          end do
        end if
      end do
      
      do idom = 1, mesh%ndom
        if (mesh%domain(idom)%active) then
          call mesh%get_params(idom=idom, ids=ids, ide=ide, jds=jds, jde=jde)
          associate(dzsdxL=>static%domain(idom)%dzsdxL, &
                    dzsdxR=>static%domain(idom)%dzsdxR, &
                    dzsdxB=>static%domain(idom)%dzsdxB, &
                    dzsdxT=>static%domain(idom)%dzsdxT, &
                    dzsdyL=>static%domain(idom)%dzsdyL, &
                    dzsdyR=>static%domain(idom)%dzsdyR, &
                    dzsdyB=>static%domain(idom)%dzsdyB, &
                    dzsdyT=>static%domain(idom)%dzsdyT  )
            !$OMP PARALLEL
            !$OMP DO COLLAPSE(2)
            do j = jds, jde
              do i = ids, ide + 1
                ! unify dzsdx on x-dir
                dzsdxL(:,i,j) = 0.5 * ( dzsdxL(:,i,j) + dzsdxR(:,i-1,j) )
                dzsdxR(:,i-1,j) = dzsdxL(:,i,j)
                
                ! unify dzsdy on x-dir
                dzsdyL(:,i,j) = 0.5 * ( dzsdyL(:,i,j) + dzsdyR(:,i-1,j) )
                dzsdyR(:,i-1,j) = dzsdyL(:,i,j)
              end do
            end do
            !$OMP END DO
            !$OMP DO COLLAPSE(2)
            do j = jds, jde + 1
              do i = ids, ide
                ! unify dzsdx on y-dir
                dzsdxB(:,i,j) = 0.5 * ( dzsdxB(:,i,j) + dzsdxT(:,i,j-1) )
                dzsdxT(:,i,j-1) = dzsdxB(:,i,j)
                
                ! unify dzsdy on y-dir
                dzsdyB(:,i,j) = 0.5 * ( dzsdyB(:,i,j) + dzsdyT(:,i,j-1) )
                dzsdyT(:,i,j-1) = dzsdyB(:,i,j)
              end do
            end do
            !$OMP END DO
            !$OMP END PARALLEL
          end associate
        end if
      end do
    end associate
  end subroutine topo_unify_edges_deriv

  subroutine topo_final()

    if (allocated(topo_lon     )) deallocate(topo_lon     )
    if (allocated(topo_lat     )) deallocate(topo_lat     )
    if (allocated(topo_zs      )) deallocate(topo_zs      )
    if (allocated(topo_landfrac)) deallocate(topo_landfrac)

  end subroutine topo_final
  
end module topo_mod
